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, We study a generalized honeycomb lattice spin- 1/2 Heisenberg model with nearest-neighbor anti- 

T-H ' ferromagnetic 2-spin exchange, and competing 4-spin interactions which serve to stabilize a staggered 

' dimer state which breaks lattice rotational symmetry. Using a combination of quantum Monte Carlo 

(N : numerics, spin wave theory, and bond operator theory, we show that this model undergoes a strong 

' first-order transition between a Neel state and a staggered dimer state upon increasing the strength 

' of the 4-spin interactions. We attribute the strong first order character of this transition to the 

spinless nature of the core of point-like Z3, vortices obtained in the staggered dimer state. Unlike 
in the case of a columnar dimer state, disordering such vortices in the staggered dimer state does 
not naturally lead to magnetic order, suggesting that, in this model, the dimer and Neel order 
parameters should be thought of as independent fields as in conventional Landau theory. 



in 



^ \ I. INTRODUCTION 

^ ■ 

, Phase transitions between two phases of matter distinguished by symmetry properties are usually well-described by 
Landau theory. In its simplest form, Landau theory expresses the free energy of the system as an analytic function 
of the order parameter whose value captures the symmetry breaking inherent in the ordered state. All analytic 
' terms consistent with the symmetries of the microscopic Hamiltonian are included in this free energy function, with 
, coefficients that are undetermined functions of the microscopic interactions in the system. The phase of the system in 
(— I this description is obtained by minimizing this free energy function over different values of the order parameter. In this 
Q description, phase transitions to the broken symmetry long-range ordered phase are driven by changes in the values of 
O . the coefficients of various terms, which change the position of the minimum to a non-zero value of the order parameter 
at the phase transition. Symmetry considerations, which dictate the form of the terms allowed in the Landau free 
. energy function, then allow one to decide whether a particular phase transition is generically a first order transition or 
J> ' second-order in nature. For instance, for systems with a global Z2 symmetry and a scalar order parameter (such as the 
Ising model), this Landau theory approach correctly predicts that the transition to the symmetry breaking long-range 
ordered phase is generically a second-order transition, with the order parameter growing continuously from zero at 
the transition. On the other hand, when the two phases on either side of the transition break different symmetries, 
and thereby possess different order parameters, the prediction of Landau theory in the generic case is that the change 
' from one phase to the other proceeds either by via multiple transitions (going through an intermediate phase with 
coexisting orders or with both orders being absent), or via a direct first-order transition with one order parameter 
' abruptly jumping to zero and the other abruptly becoming non-zero at precisely the same transition point. 
. . , Recent work by Senthil and co-authorsi*^ suggests that such a Landau theory approach is misleading for a class 
^ • of quantum phase transitions between Neel ordered antiferromagnets and valence-bond ordered paramagnets, most 
Jv>( , notably in two-dimensional square-lattice antiferromagnets. In this case. Landau theory would proceed by writing 
j_j ■ down the free energy function as an expansion in the Neel order parameter n and the valence-bond order parameter if). 
' Since one of these lives in spin-space, and the other represents order in real-space. Landau theory considerations would 
predict a first-order transition or an intermediate phase in the generic case. Such an intermediate phase appears likely 
for instance on the honeycomb lattice, when an additional next-nearest neighbour exchange coupling destroys the Neel 
ordering of the nearest neighbour Heisenberg antiferromagnets. However, RefQ argues that the transition can in fact 
be a generically continuous transition, and is better described in terms of 'deconfined' spinon variables rather than the 
order parameter fields of Landau theory. This failure of Landau theory has been ascribed to the presence of crucial 
Berry phase terms in the action written in terms of the order parameter fields, and it is conjectured that when these 
are correctly taken into account, one arrives naturally at a description in terms of 'deconfined' spinous interacting 
with a [/(I) gauge field. This continuum description leads to a prediction of a direct second order transition between 
the two phases if certain monopole operators allowed in the theory are actually irrelevant at the fixed point describing 
the transition, rendering the emergent gauge field effectively non-compact. This was conjectured to be the case for the 
transition from the Neel ordered antiferromagnet to a four-fold symmetry breaking valence bond solid phase in square 
lattice antiferromagnets. Various numerical works on a particular spin model with multiple-spin interactions appears 
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to lend support to these field theoretic ideas. Furthermore, even in cases where these monopole operators are not 
irrelevant, this approach suggests that the transition would be weakly-first order, allowing one to use this continuum 
description to describe the physics at all but the largest length scales at which the weakly-first order nature of the 
transition asserts itself. 

A more intuitive view of 'deconfined criticality' was provided in subsequent work by Levin and Senthil^ who came up 
with a simple picture for the spin- 1/2 spinon variables and the C/(l) gauge field that make up the basic ingredients for 
the correct description of a generically continuous transition between a Neel ordered antiferromagnet and a columnar 
valence-bond solid (VBS) ordered paramagnet on the square lattice. The basic idea was that a Z4 vortex formed by 
breaking up the sample into four domains of solid order meeting at the vortex core necessarily contains a free spin-1/2 
variable localized at its core. The core energy of these vortices is expected to decrease upon approaching the vicinity 
of the transition to the antiferromagnetic phase. Under the assumption that the anisotropy is irrelevant, the phase 
of the order parameter winds continuously around the core (as in U{1) vortices). It is then natural to think of this 
transition to the antiferromagnet in terms of the proliferation of these C/(l) vortices, and write down a theory for the 
transition in these vortex variables. Since solid order is destroyed by the proliferation of these vortices, the destruction 
of the solid order is accompanied by establishment of spin ordering corresponding to the condensation of these spin-1/2 
degrees of freedom in the cores of the proliferating vortices. Reasoning in this manner. Levin and Senthil were able 
to deduce the form of the continuum theory for the resulting phase transition from symmetry arguments and these 
simple intuitive considerations. 

Motivated by this simple argument, it seems reasonable to make the following supposition: If the lowest core-energy 
vortices in the valence-bond order have a free spin in their cores, then the transition to the adjacent antiferromagnetic 
phase may be expected to admit a natural description in spinon variables, of the type developed in Ref. [l]. This would 
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FIG. 1: A four-spin operator in the JQ model on the honeycomb lattice (left panel) and a six-spin operator (right panel)in the 
staggered JQ3 model on the square lattice. 



imply either a generically continuous phase transition if the Zn spatial anisotropy (with n = 4 for the square lattice) 
is irrelevant at the critical fixed point, or a weakly-first order transition if this anisotropy turn out to be relevant. On 
the other hand, if the lowest core-energy vortices have no free spin in their core, there is no natural way to obtain 
magnetic order by proliferating vortices in the dimer order; one may then expect standard Landau theory to be valid, 
and the transition to a nearby antiferromagnetic state would then proceed via an intermediate phase or be strongly 
first order. 

Some evidence for this point of view has already been provided by recent work on models that exhibit a strongly 
first order transition between a 'staggered' valence bond solid order and Neel order on the square lattice^i^ in a model 
with nearest neighbour exchange and competing 6 spin interactions on the square lattice: 

H ^~jY,P,j-Q:i P^PkiPmn, (1) 

(ij) {ij,kl,"in) 
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FIG. 2: Cartoon of vortices in staggered valence bond nematic on the honeycomb lattice (left panel) and the staggered valence 
bond solid (right panel) on the square lattice. Note that one is not forced to have a free spin in the core of these vortices. 



where denotes a bipartite singlet projector, Pij = -j — • S^. The plaquette interactions Q3 in the formula above 
are represented pictorially in Fig. [1] (right panel) . 

As is clear from Fig [5] (right panel), the simplest caricatures of vortices in the staggered valence bond solid order on 
the square lattice are indeed without any free spins at their core. In the staggered case, it is also possible to construct 
vortices with spins at their cores, but these are expected to cost more energy due to additional singlet formation 
cost associated with leaving one spin free in the core, as is clear from Fig |3] (right panel) . This suggests that the 
transition would be strongly first order, which is consistent with recent numerical works^. This should be contrasted 
with studies of very similar "JQ" models on the square lattice, with the plaquette interactions chosen to favour 
columnar order—. Extensive numerical work on these models has led to the conclusion that the transition from Neel to 
columnar valence bond solid order is indeed a continuous transition^^^— although the critical point exhibits logarithmic 
violations of standard finite size scaling^iii. These violations of scaling have also been interpreted as possible evidence 
of first order behaviour . SU{?>) and SU{A) versions of the same transition have also been studied^ii^ii^ and found 
to be continuous in nature. 

Here we focus on a related system on the honeycomb lattice: On the honeycomb lattice, the staggered ordering 
of near-neighbour valence bonds does not break lattice translation symmetry, but does break symmetry of three-fold 
lattice rotations. As is clear from Fig [2] (left panel), the simplest construction of elementary vortices in this 
valence-bond nematic order naturally yields vortices with spinless cores. Again, vortices with spinful cores are also 
possible, but are expected to cost more core-energy due to two reasons. First, forming a spinful vortex core shown 
in Fig [3] (left panel) involves the breaking of additional singlet bonds, and would cost an additional core energy of 
the order of the spin gap. At the same time, the spinless vortex core, shown in Fig. [2] (left panel), allows for dimer 
resonances on the core plaquette which is expected to lead to core energy lowering. 

The transition between such a valence-bond nematic and a nearby antiferromagnetic phase on the honeycomb lattice 
is thus expected to be strongly first order based on our supposition above. Here we focus on a spin Hamiltonian with 
nearest neighbour exchange couplings that compete with four-spin interaction terms that stabilize a valence-bond 
nematic state. We study this system numerically by a projector Monte Carlo algorithm in the valence-bond basis, as 
well as analytically by spinwave and bond-operator expansions. Our numerics finds a strongly first order transition 
between a valence-bond nematic and a antiferromagnet, consistent with the intuitive picture outlined above. In 
addition, we use our approximate analytical calculations to provide a reasonable semi-quantitative account of our 
numerical results. We note that an exactly soluble model with such a valence bond nematic ground state has also 
been proposed recently*^ 

The remainder of the paper is organized as follows: In section |TT1 we define the model studied, and describe the 
approximate analytical methods that we use for their study. In section lllli we describe our numerical approach and 
summarize the results obtained. Finally, we conclude with a brief discussion in section ITVl 



FIG. 3: Vortices with free spins in their cores are also possible in staggered valence bond nematic on the honeycomb lattice 
(left panel) and the staggered valence bond solid (right panel) on the square lattice, but are expected to cost higher core energy 
than ones with spinless cores. 



II. MODEL AND ANALYTICAL TREATMENT 



The honeycomb lattice staggered JQ model Hamiltonian is defined as folows, 

H=-JY,Pv-Q J2 P-oPki , (2) 

{ij) {■ijM) 

where the four spin plaquette operators are products of two bipartite singlet projectors acting on two parallel bonds 
of the hexagons in the honeycomb lattice (see Fig.[ljleft panel)). This Hamiltonian can be written explicitly as a sum 
of two-spin and four-spin interactions as follows, 

HjQ = H2 + Hi (3) 
H2 = J^S, -S, (4) 



Ha = - 



QY, [(S? . S?)(S? . S«) + (S? . s«)(sp . S«) + (S? . s?)(sp . S?)] (5) 

o 

where denotes a sum over all hexagons of the lattice, SP denotes the i-th spin on a hexagon, with i — 1-6 
labelling the 6 sites clockwise around the hexagon, and J = J + Q/2. 

At a heuristic level, we expect that the Q term drives a transition to a valence bond nematic state with staggered 
ordering of near-neighbour singlet bonds, while the J term favours a Neel state with spins on the A-sublattice [B- 
sublattice) of the honeycomb lattice aligned (anti-aligned) along a spontaneously chosen axis in spin space. A simple 
approximate description of the Neel phase is obtained by doing a standard spin wave analysis. On the other hand, to 
describe the valence-bond nematic, we use a description in terms of bond-operators that becomes very accurate in the 
limit of strong nematic order. By comparing the predictions for the ground state energy from these two approximate 
calculations, we also estimate the value of Q at which we expect a transition between the two phases. 

A. Spin wave theory of the Neel state 

In the Neel ordered phase of the JQ model defined above, let us assume we can decouple the four spin interactions 

as 

(S? . S?)(S? . SP) = -a(S? • + S? . S?) - a" (6) 
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where a = — (S; • Sj) is the nearest neighbor spin correlation which is assumed to be the same on all neighboring 
bonds and which needs to be determined self-consistently. Such a decomposition can alternatively be obtained within 
a path integral formalism by using a Hubbard Stratonovich (HS) decoupling of the quartic spin operators followed by 
a static mean field theory of the HS field. In this case, the Hamiltonian reduces to the form 



^mf = Jeff J2 ^» • + SQa^iVn (7) 

where we have defined Jeff ~ J + 2Qa, and Nh is the number of unit cells on the lattice (the number of spins is 
Ns = 2Nh). 

We treat this effective Hamiltonian using the standard spin wave theory by defining Holstcin Primakoff bosons via 
{S\, S^) = (i — a^a ,a ,a'') and {Sg, S^, S]^) = (6^6 — ^,b\b ) for sites on the A-sublattice and B-sublattice 
respectively. Retaining terms to quadratic order in the bosons leads to the Hamiltonian 

where Fk = (1 + e"*'^'' + e"'*^""'^'') with fc^ = k • a and = k ■ b. Setting Fk = iFkle*^'', we can diagonalize this 
Hamiltonian via a Bogoliubov transformation 



/cosh6'ke*'^'' -sinh^k 
blj " ^ -sinh6lk cosh^ke"''^''; \d_^^ 

where sinh2(?k = iTkl/i^k, cosh 20k = S/ilk, and fik = ^9 — |Fkp. This yields the diagonalized Hamiltonian 



(9) 



i/mf = -JjeffA^// + 3Qa27V^ + :^E(r!k-3) + :^Er!k(4ck + 4^k) (10) 

k k 

with a ground state energy 

k 

Applying the Feynman-Hellman theorem to the mean field Hamiltonian, _ffmt in EqlTl yields 

1 dEsw 1 1 V- 



3Nh 9Jeff 4 6Nh V 

k 



(r!k-3) (12) 



Eqns. (fTTj) and (fT2|) . together, determine the ground state energy of this model in the Neel ordered phase. A feature 
of our mean-field spin-wave analysis is that a is independent of Q/ J. Numerically, we find a « 0.3549 within this 
spin wave approach which agrees with earlier spin wave results^®. This value is quite close to the value, a sa 0.3627, 
deduced from a recent optimized valence bond trial wave function studj^ of the nearest neighbor Heisenberg model, 
and series expansion results^* on the Heisenberg model which give a « 0.3659. This spin wave result for a must be 
used as input to compute the ground state energy per spin of the J-Q model which is given by egw = — |(J + Qa)a; 
this is plotted in Fig|4l Within this approach, the ordered moment is also independent of Q/ J, and is given by 

k 

and takes on a value m ~ 0.2420, whereas series expansion studies of the nearest neighbor Heisenberg modeti^ yield 
m 0.266. In the next section, we will check to what extent the numerical results bear out our prediction of a nearly 
Q-independent value of m in the Neel phase. 



B. Bond operator theory of the staggered dimer state 

In the other limit, when the Q term dominates over J Hamiltonian at large Q has spontaneous staggered dimer 
order in the ground state, we can use a bond operator formalism to compute the ground state energy and correlations 
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in this phase. We label the spins by their unit cell position r = ma + nh and a sublatticc index p = 1,2. In the 
staggered dimer state, let us assume that the spins at sites (r, 1) and (r, 2) in every unit cell form a singlet. This dimer 
pattern then spontaneously breaks the lattice rotational symmetry but leaves the translational symmetry intact. We 
define the singlet and triplet states on this pair of sites in terms of bond operators via 

(14) 
(15) 
(16) 
(17) 

where we must satisfy the constraint s\sj. + 4c«*ra = 1 at each site. The spin operators can be rewritten in terms of 
the bond operators as 

= li-^r [4tra + tLSr] ' '^e^P^^t^j (18) 

where a = x,y,z. Assuming the staggered dimer state corresponds to a uniform singlet condensate, we can replace 
sj — > s and — > ,s. Furthermore, if the singlet condensate is robust, so that the dimer order is strong as observed in 
the QMC numerics, we expect the triplet density to be small and a Bogoliubov-type theory, where the triplet density 
plays the role of a small parameter, to be a reasonable starting point. Expanding out the Hamiltonian in terms of 
these bond operators, and retaining leading quadratic terms in the triplet operators, we find 



\s) = 




- Ilt2) 




\z} = 




+ ilt2) 




\x) = 




- 4-14-2) 




\y) = 




+ 4.42) 
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- J-ryitla+traW t +t f - i 

A £. — icK/\ r— o.a v—h.OL r—a—b,a r 



4 ^ ^ ra/v r— tijO; r—b.a r—a—b^ot r—a—b,ct 

r 

- Ps^Nh--^Qs^Nh-XNh{s^-1) (19) 

where A is a chemical potential which serves to satisfy the constraint equation on average by fixing (4a*rc«) = 1 — s^- 
Going to momentum space, the Hamiltonian takes the form 



(20) 



where 



= -i^-s' + ^-^s^)Nh - XNnis' - 1) (21) 
4 16 

A = (^ + - A) (22) 

Bk = -^s^(cos kb + cos{ka + kb)) (23) 

and the prime on the momentum sum indicates that we only sum over half the Brillouin zone (keeping states with 
fco > for example). Diagonalizing this via a Bogoliubov transform 



*ka \ ^ f cosh6'k -sinh0k\ / 
t-U l-sinh^k cosh^k ; UU^ 

we find that the diagonal Hamiltonian takes the form 



(24) 



H = + E^k(cLck„ + 4.4„ + 3) - ^N^d + ^«-^) (25) 
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Q/J 

FIG. 4: (Color online) Energy per spin (in units of J) of the Neel state obtained in spin wave theory (SW, red-solid line) 
compared with the staggered dimer ordered state obtained using bond operator theory (BO, blue-dotted line). There is a 
strong first order transition at Q/J ~ 1.35 (Q/J ~ 4.2). 



where = v^A(j4 + 23^). This leads to a ground state energy 

E^o = E('^ +3^2^^- In^C^ + ^fs'). (26) 

k 

This energy must be minimized with respect to s once A is determined by the number equation 

3 ^(A±^_i) = i_,. (27) 



k 



s 
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which follows from demanding dE-Qo/dX — in order to satisfy the constraint relation (tl^t^^) — 1 

Solving these equations numerically, we self-consistently determine s and A. Using these, we compute the ground 
state energy and excitation spectrum of the valence bond nematic state. As seen from Fig. |4l we find that the energy 
of this state lies below the antiferromagnetically ordered state for large Q 7 as we expect, but there is an energy level 
crossing at Q/J ~ 1.35 (Q/J ~ 4.2), so that the Neel state has lower energy at smaller Q. (Within the bond operator 
approach, we find that the triplons condense around Q/ J ^ 0.2 which is well below the first order transition point.) 
Close to the transition, in the spin gapped phase, we find the spin-spin correlation on the dimerized bond to be 
~ —0.73 which corresponds to nearly complete dimerization, so that this is a very strong first order transition. We 
next turn to a numerically exact quantum Monte Carlo study of this model. 



III. NUMERICAL STUDY 



We use the valence bond projector Monte Carlo technique^^ to study various ground state properties of the staggered 
version JQ model on honeycomb lattice as a function of coupling Q/J. Taking advantage of the improved loop- update 
technique developed in Ref. [lO, we scan the phase diagram of this model for systems with upto 2 x 32 x 32 sites. A 
Monte Carlo projection length of 6L^ is used to ensure the convergence of observables to the ground state expectation 
values (we have checked the convergence by comparing with the results of exact diagonalization studies at small sizes) . 
Also, to counter ergodicity problems in the Monte Carlo simulations in the valence bond nematic phase and in the 
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vicinity of the first order transition, we employ many different Monte Carlo moves to update the bond and plaquette 
operators in a Monte Carlo configuration. These updates include attempts to reflect a plaquette operator about its 
diagonal bond (FiglSK), to rotate the diagonal plaquette operators about a site (FiglSlD), to transform a plaquette 
operator into a bond operator by deleting a bond and vice versa (Fig|SJ:), and so on. 





C) 

^_ (1/4-Si^S^) (1/4-Si^S^) or -(S.. sM 

FIG. 5: Monte Carlo moves to update plaquette operators: a) reflecting a plaquette operator about a bond with a diagonal 
operator acing on it, b) rotating a diagonal plaquette operator about a site, and c) transforming a diagonal or off-diagonal 
bond operator into a plaquette operator and vice versa. 

We find it convenient to define the valence bond nematic order parameter, ■ip{r) = S(r) • S(r + ei) +exp (i27r/3)S(r) • 
S(r + 62) + exp (i47r/3)S(r) • S(r + 63), where ei 2,3 denote the three nearest neighbor bond vectors of the honeycomb 
lattice. By monitoring the correlations of "0 simultaneously with those of the spins, we are able to distinguish easily 
between the Neel antiferromagnet, and the rotation-symmetry breaking valence bond nematic. To do this we measure 



^(1) = 2i^E("W 1))' 



(28) 
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where n(r) = S(r) — S(r + ei), and the sum over r runs over all A-sublattice sites of the lattice. We also define 

m^^J2(^*irWr + l)), (29) 

r 

where r = (x, y), and 1 — {L/2, L/2) with L being the linear dimension of the system and Ng ~ 2L^ being the number 
of sites in the system. In the presence of long range antiferromagnetic order, the asymptotic large size limit of C(l) 
tends to a finite value, which equals the square of the staggered magnetisation 

M' = ^ E (30) 

r£A,r'£A 

In the absence of long range antiferromagnetic order, we expect this quantity to vanish in the large size limit. Similarly, 
D{1) serves as a order parameter for long range valence bond nematic order. 



A. Results for honeycomb lattice model 




2 4 6 8 10 

Q/2J 



FIG. 6: Simultaneous first order jump in Neel order parameters C(l) and valence bond nematic order parameter D{1), as Q/ J 
is varied for the JQ model on the honeycomb lattice. 

From our numerical study of these observables, we deduce that a transition from the Neel phase to the valence bond 
nematic phase takes place as Q/J is increased to the vicinity oi Q/J ^ 7. In this vicinity, there are clear first order 
jumps in both the order parameters as shown in FiglHl At first sight, the data suggests that the transition point, as 
determined from the location of the jump, may have a tendency to drift with increasing system size. However, this is 
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actually an artifact of the slowness of the code near the first order transition region at larger sizes, as is clear from an 
alternate more precise analysis which we now summarize. In this alternative approach, we determine the transition 
point by locating the kink in the plot of energy as function of tuning parameter at smaller sizes, where there are no 
such issue of slowness of the Monte Carlo update procedure used for calculating ground state properties. As shown 
in Fig. [71 the energy per site shows a linear behaviour, but with slopes of the line changing discontinuously near 
the transition. This can be understood as an avoided level crossing taking place at the first order transition point. 
Identifying the crossing point of these two linear extrapolations for the energy on is a convenient way to determine the 
transition point^ with reasonable accuracy. Using this method, with data from three different system sizes L = 12, 16 
and 24, we identify the first order transition to be at Q/ J — 6.4 ± 0.2. 

0.015 
0.01 
0.005 

o 

^ 

+ 

-0.005 

P -0.01 
c 

^ -0.015 
-0.02 
-0.025 

4.2 5.2 5.6 6 6.4 6.8 7.2 7.6 8 

Q/J 




FIG. 7: An avoided level crossing of the two competing candidate ground states, showing up as a kink in the ground state 
energy of the honeycomb lattice JQ model, when plotted as a function of the tuning parameter Q. This serves as clear signature 
of the first order transition taking place in this model system. Note that the energy difference from the fitted linear behaviour 
on the antiferromagnetic side is plotted, not the absolute ground state energies on both sides - this is just to highlight the 
discontinuity of slopes of the two curves at the transition point. To obtain linear fits to energy values on the valence bond 
nematic side, the data points at Q/J < 6.8 have not been considered. 



First order transitions are associated with two coexisting free energy minima corresponding to the two phases. 
This leads to the appearance of hysteresis effects in the vicinity of the transition for all but the most efficient Monte 
Carlo algorithms when a configuration from one of the phases is evolved by tuning the coupling constant across the 
critical value at finite rate (i.e without allowing for an extremely large number of equilibriation steps between each 
successive change of coupling constant). For example, starting with a configuration deep in the valence bond nematic 
phase, the system remains stuck in the metastable valence bond nematic state even when Q is reduced to less than 
Qc and conversely, the system remains in the Neel phase upon ramping Q up beyond Qc if one starts with a typical 
configuration obtained from the projection algorithm deep in the Neel phase. We present such a hysteresis plot for 
antiferromagnetic order parameter (FigjS]) near the transition to emphasize the first order nature of the transition. 

Our QMC simulations also yield a characteristic double peaked order parameter histogram near Qc as shown in 
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Q/2J 

FIG. 8: Hysteresis of (M^) around transition point for honeycomb lattice JQ model. 

Fig. [S] However, the 'time' (number of Monte Carlo steps) taken for the system to 'tunnel' between the two phases 
is extremely large, and it is therefore prohibitively expensive to perform runs that could encompass a large number 
of such tunneling events. As a result, it is not possible with our limited computational resources to obtain accurate 
estimates of the relative heights of the two peaks in the histogram. 

Finally, we note that in a variant of staggered JQ2 model on honeycomb lattice, where parallel bonds from neigh- 
bouring hexagons make up the plaquette interactions, we find that increasing Q fails to destroy AF order. For this 
model, as shown in Fig. [TUl there is no hint of a transition even at large Qj J . 

IV. DISCUSSION 

We have thus studied an example of a transition between a Neel ordered antiferromagnet and a three-fold rotation 
symmetry breaking valence bond nematic on the honeycomb lattice. Our numerical results showed that the transition 
is strongly first order, consistent with the general scenario outlined in the introduction. The transition point obtained 
from our analytical considerations is at Q/ J w 1.35 [Q/ J ~ 4.2), which is in reasonable agreement with the results 
of our QMC simulations, which give a strong first order transition at Q/J « 1.5 {Q/J « 6.4). Furthermore, the 
nearly Q independent values of the staggered magnetization Ad in the Neel phase, and of the nematic order \D\ in 
the staggered dimer state, are also consistent with our analytic considerations. These observations suggest that the 
bond operator approach provides a reasonable account of the dimer phase in situations where the kinetic fluctuations 
of dimers are not important in the vicinity of the transition. This appears to be different from the Neel to columnar 
dimer transition in the square lattice case where such singlet fluctuations appear to be important in driving the dimer 
order to zero at the quantum phase transitioni^i 
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L=24, Q/2J=3.75 
Q/2J=4.00 
Q/2J=3.00 
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FIG. 9: Histogram of valence-bond-nematic order parameter \D\ = \D\{1) for L=24 at three different Q values for honeycomb 
lattice JQ model. Presence of double peak structure at Q/ J = 7.5 indicates first order nature of transition. 
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